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Abstract 

We present an automated generation of the subtraction terms for next-to- 
leading order QCD calculations in the Catani-Seymour dipole formalism. For 
a given scattering process with n external particles our Mathematica package 
generates all dipole terms, allowing for both massless and massive dipoles. The 
numerical evaluation of the subtraction terms proceeds with MadGraph, which 
provides Fortran code for the necessary scattering amplitudes. Checks of the 
numerical stability are discussed. 
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Numerical evaluation of subtracted matrix elements interfaced to MadGraph (3]-[5]| (stand-alone 
version) using helicity amplitudes and the HELAS library (contained in MadGraph). 

Restrictions on complexity of the problem: Limitations of MadGraph are inherited. 

Typical running time: Dependent on the complexity of the problem with typical run times of the 
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1 Introduction 

The Large Hadron Collider (LHC) allows us to explore an energy regime far beyond what has 
been accessible in direct measurements up to now. Operating at the TeV scale it will shed light 
on the mechanism of electroweak symmetry breaking. The large center-of-mass energy allows the 
production of new heavy particles — if they exist. It is expected that the results obtained from the 
LHC will significantly influence our future understanding of nature. The high potential of the LHC 
in making an important step forward towards a deeper understanding of elementary particle physics 
comes not for free. Experimentally as well theoretically the LHC is a very challenging experiment. 
The large multiplicity in the individual event together with pile-up and the complications due to the 
underlying event make the experimental analysis at the LHC highly non-trivial. The separation of 
known Standard Model (SM) physics from new physics is very demanding, in particular since in 
many cases the new physics signals are overwhelmed by the large SM backgrounds. Sophisticated 
methods have been devised in the past to cope with this situation within the experiments with 
respect to theoretical predictions. 

Colliding protons at LHC interact primarily through strong interactions and Quantum Chromo- 
dynamics (QCD) plays an important role in this context. It is well known that perturbative QCD 
predictions are in general plagued by a large (residual) renormalization and factorization scale de- 
pendence. It is not rare that predictions for cross sections change by 100% when the scales are 
varied in a reasonable range. To reduce the large scale uncertainty next-to-leading order (NLO) 
calculations are required, which generally consist of two ingredients: one contribution from the 
virtual corrections and a second contribution from real emission of one additional parton. With 
respect to the virtual contributions much progress has been achieved recently, although the evalu- 
ation of one-loop amplitudes for "large" multiplicities (2 — » 3, 2 — » 4) is still a highly non-trivial 
enterprise. Fortunately, as far as the real corrections are concerned the situation is much better 
and efficient methods exist for the numerical evaluation of the required matrix elements. However, 
the integral over the phase space leads to soft and collinear singularities, which we denote here 
generically as infrared (IR) divergencies. In combination with the virtual corrections, all IR cancel 
between the two contributions for the physical observables of interest [|8l-fT0l. 
General algorithms, typically classified either as slicing or subtraction methods, are available for 
the extraction of soft and collinear singularities encountered in the real corrections. In both cases 
one makes use of the universal behavior of QCD amplitudes for soft and collinear configurations. 
IR singularities are nowadays usually regulated within dimensional regularization, so the same 
regulator has to be applied to the real corrections in a way that numerical integration over the 
phase space in 4 dimensions is possible in the end. 

In slicing methods the idea is to separate the phase space into resolved and unresolved contribu- 
tions [fTTT - Tl4l . Unresolved regions at NLO are those where one parton becomes soft or two become 
collinear. In these regions the matrix elements are approximated using the QCD factorization the- 
orems. After this simplification the unresolved regions can be integrated analytically in d = 4 - 2e 
dimensions and the emerging singularities cancel analytically against the corresponding ones in the 
virtual corrections. Resolved regions on the other hand exclude all IR singularities by definition 
and can be integrated numerically in 4 dimensions. A one-dimensional illustration of the slicing 
approach is shown below, 
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= i/(0) + ln(<5)/(0)+ f ^Jx + 0(6), (1.1) 

where /(jc) is an arbitrary function, which is regular at x = 0. It is a drawback of the slicing 
method that partial results exhibit a logarithmic dependence on the cut 8 which separates resolved 
from unresolved configurations. This logarithm — which cancels when the two parts (resolved 
and unresolved) are combined — is manifest in the analytic integration of the unresolved terms. 
However, for the resolved terms the logarithmic dependence arises from the numerical phase space 
integration. Since the matrix elements in the unresolved phase space regions are only approximate 
one tends to make the respective region around singular configurations as small as possible. This 
procedure, though, would result in large numerical cancellations and a loss of accuracy in the sum 
of the resolved and unresolved parts compared to the accuracy reached in the numerical integration. 
In practice this requires a compromise between the quality of the approximation and the numerical 
effort to achieve a certain precision in the sum of both, resolved and unresolved contributions. 
Subtraction methods make use of our knowledge about QCD factorization in the soft and collinear 
limits to construct suitable "counter- terms". These have to match pointwise all singularities in the 
real-emission matrix elements and, at the same time, should be simple enough to be integrated 
analytically in d = 4 - 2e dimensions over the entire phase space 03 13 [21 • A one-dimensional 
example is shown in Eq. (11.2b . 



Here the cancellation of the singularity takes place at the integrand level — i.e. the divergencies 
cancel pointwise — and the integrals are easier to evaluate numerically. However, care must be 
taken with respect to the numerical accuracy. Deep in the singular regions the individual con- 
tributions become arbitrary large and, due to limited numerical precision of floating point arith- 
metic, their cancellation might be incomplete leading to potentially wrong results. While our 
one-dimensional example in Eq. (1 1 .21) admits the construction of a suitable subtraction term in a 
straightforward manner this no longer true when considering complicated scattering amplitudes. 
Fortunately, this formidable problem has been solved with the Catani-Seymour dipole formal- 
ism 03I21- (Similar algorithms are presented in Refs. 03H21-) Based on the factorization of soft 
and collinear singularities in an SU(A0 gauge theory all subtraction terms are constructed from 
universal functions and process specific amplitudes. It turns out, though, that for complicated 
processes the complete subtraction term is a sum over many different contributions — tedious to 
derive by hand. On the other hand, the method is completely algorithmic, thus an automation is 
feasible. This is a timely problem and its solution is the aim of the present work IfTTlfTHl . We note 
that it has recently also been addressed by other groups [[T9l - l2~2"l . see also 11231 . 
In passing let us briefly mention that the attractive features of the subtraction approach (universality 
of counter-terms and numerical stability) have led to extensions of the formalism to next-to-next- 
to-leading order (NNLO), see e.g. Il24l - l27l . Presently, the proposed schemes at NNLO apply to 
processes without colored partons in the initial state and an arbitrary number of massless parti- 
cles (colored or colorless) in the final state. The necessary counter-terms are either derived from 
so-called antenna functions ll28ll or alternatively defined as universal counter-terms based on QCD 
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factorization in the various soft and collinear limits of singly- and double-unresolved parton con- 
figurations ||29] - l32l . Ongoing work at NNLO is concerned with extensions to colored partons in 
the initial state as needed for the LHC. 

The outline of the paper is as follows. In Section [2] we review the general features of the Catani- 
Seymour algorithm HI 13 and in Section [3] we describe the details of its implementation in the 
AutoDipole package. Particular emphasis is put on the details of the numerical evaluation of sub- 
traction terms via an interface to MadGraph (3H51, which uses the HELAS library (6113 for the 
computation of helicity amplitudes. Section |4] illustrates with a few examples of how to use Au- 
toDipole in practice. The examples serve also as a non-trivial cross check of the implementation 
when compared with existing results from the literature. Finally, we conclude in Section [5] and 
list some technical details of the implementation and the comparison with the literature in Appen- 
dices IaHS 



2 Review of the Catani-Seymour subtraction formalism 

We consider a generic scattering process for the production of an n-parton final state. At NLO 
accuracy in QCD the corresponding cross section may be written as: 

o"nlo = o"lo + £cr NLO , (2.1) 

where cr LO denotes the Born contribution at leading order (LO) and the genuine NLO correction 
£cr NLO receives contributions from three different sources, 

8cr NLO = I dcr vin + I Jcr real + I dx I dcr iact . (2.2) 

Jn Jn+l J Jn 

By dcr vin we denote here the contributions from the virtual corrections and by <icr real the ones from 
the real emission of one additional parton. For hadrons in the initial state there is a third contribu- 
tion rfcr fact due to the factorization of initial state singularities. The subscript on the integral signs 
in Eq. (|2.2I) indicates the dimensionality of the phase space: The virtual corrections are integrated 
over an n-parton phase space while the real emissions are integrated over an (n + l)-parton phase 
space. All three contributions are individually divergent and thus require regularization in inter- 
mediate steps until the divergencies are canceled. In the Catani-Seymour formalism the expression 



for £cr NLO in Eq. (12.21) is rewritten schematically: 

8<t nlo = I (<icr real + dA) + I (dcr virt + I dA')+ I dx I (Jcr fact + dA"), (2.3) 

Jn+l Jn J I J Jn 

with 

0=f dA+ f fdA' + \ dx fdA" . (2.4) 

Jn+l JnJl J Jn 

The expressions dA, dA' and dA" are defined such that they render the individual pieces in Eq. (12.31) 
finite. In the Catani-Seymour formalism, they correspond to the sum of all dipoles (dA), the 
integrated dipoles (Jj dA'), I-term, and the terms arising from mass factorization (dA") including 
the so-called P-, K- and H-terms. 

The explicit form of dA is constructed from the knowledge about the soft and collinear factoriza- 
tion of QCD amplitudes, which exhibit a simple factorization for collinear configurations. The 
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expression for dA in Eq. (12.31) is obtained as a sum over potentially collinear partons since the 
factorization for soft configurations can be derived from the collinear behavior of the amplitudes. 
For soft singularities however only color-ordered amplitudes exhibit a simple factorization, which 
implies in general non-trivial color correlations at the level of the squared matrix elements. These 
color correlations are reflected by the reference to an additional (so-called) spectator parton. The 
subtraction dA is thus written as a sum of individual dipoles in the following form: 

dA = Y j D(i,j;k), (2.5) 

where the sum runs over all colored partons in scattering process and the possible configurations 
for {i,j;k} are determined from the real corrections for which the subtraction term is constructed. 
The generic form of the dipoles is given by 

D(i,j;k) = dij x (l,..Aj,...X---,n\Wij i k\l,...ij,...,k,---,n). (2.6) 

The singular behavior is contained in the pre-factor dij which is essentially the intermediate prop- 
agator before the splitting into i + j. The bra «... |)and ket (| . . . » notation for the amplitude is used 
since it appears as a vector in color and spin space. Accordingly \ij,k acts as an operator in color 
and spin space and introduces the non-trivial color and spin correlations mentioned above. 
As a concrete example we display the expression for D(gi,gj,k) with all three partons in the final 
state, i.e. 

1 ~ _ Tjfc-Ty 

D(gi,gj,k) = (1, . ..ij,...,k,.. .,n\ Vg.g. ,fc|l,. ..ij, ...,k,...,n) . (2.7) 

iPi-Pj T.j 

D(gi,gj,k) describes the situation of a collinear splitting of the gluons i and j (emitter pair) in the 
presence of a spectator parton k. Here, the T ; are the color charge operators depending on the 
partons being in the fundamental or adjoint representation of the color SU(3) (for details we refer 
to flUEll). The function V gigj! k is given by 



V" v , = (M\V g - g -k\v) 
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(2.8) 



where Zi and are some functions of the Lorentz scalars s^, and Sji and a s is the strong 
coupling. The quantity, 

(\,...ij,...,k,.. .,n\T k ■ Tjj\\, . ..ij,...,k,.. .,n), (2.9) 

is called the color linked Born amplitude squared (CLBS). The amplitude (| . . .)) is deduced from 
the real emission by the factorization of one splitting. The color factor is extended by two color 
operator insertion, • Ty. The CLBS has the reduced kinematics as input. In the example the 
original momenta (pi,Pj,pt) are reduced to (pij,pk) where for example, the pij is defined as 

lJ 1 3 l ~yij,k k 
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Generally, (n + 3) original momenta are reduced to the in + 2) momenta in each dipole term. One 
distinguished feature of the dipole subtraction is that the reduced kinematics satisfy the on-shell 
conditions and momentum conservation, which makes it possible to evaluate the reduced Born 
amplitude by existing codes for LO calculations. In the gluon emitter case like the example above, 
the second term in the square bracket in Eq. (|2.8I) introduces spin correlations of the gluon. In that 
case we need also to evaluate the CLBS where the helicity of the emitter gluon in the amplitude 
is different from the one in the conjugate. In the next section we explain how the three main 
ingredients, the generation of the dipole terms, the CLBS, and the dipoles with a gluon emitter, are 
implemented in our package. 

The other two functions Jj dA' and dA" in Eq. (|2.3I) are easily explained. The first one ( L dA') 
denotes the integrated dipoles which cancel the corresponding IR divergencies in the virtual con- 
tributions at NLO. They can be constructed from universal insertion operators I, to be sandwiched 
between the bra «...|)and ket (|...)) amplitudes of the corresponding Born process, i.e. in the 
reduced kinematics with (n + 2) momenta. The expressions for the I-terms depend on the parton 
type, e.g. for massless partons we have 

I«PU) = -T^X^^Y^r^^] , (2.11) 
2/r T(l - e) Y Tf. j£ J \2pi ■ Pj) 

where the sum runs over all parton momenta {p} in the Born kinematics (corresponding also to the 
NLO virtual corrections). The poles in dimensional regularization are contained in the function *V,-, 
that is 'Viie) ~ l/e 2 and ~ 1/6 for massless partons H). For massive partons 0, the corresponding 
function < y i (e,m J -,my) contains single poles r Vi{e,mi,m.j) ~ l/e from soft gluons and the *V ; depend 
logarithmically on the parton masses m, which screen the collinear singularity, see also lfT6l[33l . 
The formulation suggested in [2] (and used in AutoDipole) allows for a smooth interpolation in the 
limit m ; — » 0. The color operator insertions, T ; • Tv, lead again to non-trivial color correlations and 
require the evaluation of the same CLBS in Eq. (12.91) with the phase space of n final state partons. 
The last term in Eq. (12.31) {dA") abbreviates the so-called P-, K- and H-terms. They are needed for 
mass factorization of initial state singularities in hadron collisions to be absorbed in renormalized 
parton distributions and, likewise, also for the final-state singularities in processes with identified 
hadrons giving rise to scale dependence of the fragmentation functions. For instance in the case of 
initial state divergencies the P-operator takes the following form 

1 u 2 
P a ' b ({p},x,n 2 F ) = ^P a \x) — VT r T,ln— ^— , (2.12) 

where fif is the factorization scale and P ab are the standard (space-like) LO splitting functions 
containing the well-known '+' -distributions in their diagonal components. Like the I-operators 
in Eq. (12.1 II) . the P-operators act on the amplitudes of the corresponding Born process, leading 
to color correlations and giving rise to the CLBS as in Eq. (12.91) . Here the dependence of P in 
Eq. (12.121) on the parton momentum fraction x leads to convolutions when integrated with the Born 
squared matrix elements or CLBS, see Eq. (12.31) , and the implementation needs a prescription for 
'+' -distributions (see e.g. Eq. (B.26) in Oil). 

Similar definitions as in Eq. (12.121) hold for the K- and H-operators. They also contain distribu- 
tions 1/(1 - x)+ that are singular in the collinear limit and parametrize the factorization scheme 
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Figure 1 : The directory structure of the AutoDipole package. 



dependence, K for factorization in the initial state and H in the final state. These remarks conclude 
our brief review of the Catani-Seymour dipole formalism. For details the reader is referred to the 
original literature Q3I21- 



3 The AutoDipole package 

The AutoDipole package constructs all necessary subtraction terms of the Catani-Seymour for- 
malism [1,2] for a given scattering process with n colored partons in a fully automatic manner. It 
can handle both massless and massive partons and it also allows for additional (non-colored) SM 
particles in the scattering process, e.g. couplings of quarks to y, Z-, -bosons and so on. 
Let us briefly sketch the details of our implementation of the dipole subtraction formalism, because 
there is a large freedom in the way how this can be performed. For instance, in Refs. Il35l[36l 
the implementation was realized in form of two independent C/C++ libraries providing all the 
necessary functions to evaluate the dipole terms. As a slight disadvantage of this approach the 
produced code is non-local and that there is some redundancy in the calculation. In the present 
work we follow a different strategy. The main idea here is to have a code generator which will 
produce an optimized fiat code. To realize this we have constructed a Mathematica program which 
acts as such a generator and its output is interfaced with MadGraph (3]-(5]|. To that end we found the 
stand-alone version of MadGraph suitable and for the numerical evaluation of all subtraction terms, 
we use helicity amplitudes provided by MadGraph (which are based on the HELAS library flMZl)- 

3.1 Structure of the code 

The complete layout of the code of the AutoDipole package is displayed in Figure [IJ The major 
part of the source code in the package (all contained in the directory lib) is written in Mathematica. 
This includes all algorithms for the generation of the dipole terms as well as the evaluation of 
the I-, P- and K-operator insertions sandwiched between the Born amplitudes in the appropriate 
kinematics. The use of a computer algebra program like Mathematica has the clear advantage 
here that all terms are accessible to symbolic manipulation. This feature is very useful for general 
studies of the IR behavior of scattering amplitudes. 
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Figure 2: The flowchart of the execution of AutoDipole. 



In using AutoDipole we are primarily interested in the generation of Fortran code for all subtraction 
terms of Eq. (12.31) . The essential parts of the flowchart are displayed in Figure[2] From a given input 
process, the Mathematica code generates all dipole terms. The important files here are dipole . £, 
reducedm. f, and the shell script for the interface with MadGraph, where dipole . f contains all 
dipole terms except for the reduced kinematics and the CLBS. The routine reducedm. f calcu- 
lates the reduced kinematics of each dipole term as defined in Eq. (12.101) . The subsequent run 
of a patched version of MadGraph (contained in the directory patch) via the interface produces 
the files, cmatrix. £ and allcolormat.inc. The routine cmatrix. f evaluates all CLBS as in 
Eq. (|2.9I) . It includes the file allcolormat . inc, which contains all extended color matrices by the 
two color operator insertions. We use the latest version of MadGraph (MG_ME_SA_V4.4.30). All 
the process specific Fortran files mentioned above are stored in a new directory under process, 
see Figured! 

In short, the evaluation of the generated Fortran code for the dipole terms proceeds along the 
following chain, see Figure [2] The routine reducedm. £ receives a given phase space point as 
input and calculates the reduced kinematics of each dipole term. Next, the routine cmatrix. £ 
(with the reduced kinematics as input) computes the CLBS of each dipole term together with the 
extended color matrices in allcolormat . inc. Finally the routine dipole . £ returns the values 
of all dipole terms having received the reduced kinematics and the CLBS as input. All this is done 
automatically when the user executes the command . /runD (see Sections |4~2l and |4~3l) to calculate 
the subtraction term for a specific phase space point. 

The algorithm of dipole term generation by the Mathematica code is shown in Section 13.21 A 
detailed account of the calculation of the CLBS is given in Section 13.31 and some aspects of the 
spin correlations for the particular case of a gluon emitter are presented in Section I3~4l The check 
of IR safety is explained in Section 13.51 An explanation of how to run AutoDipole is deferred to 
Section |H 

3.2 Algorithm of dipole generation 

Let us start with a short description of the Mathematica code generating all dipole terms. This 
procedure makes use of the following three steps: 



8 



Dipole 1 



(3) 5, 6 
(4) 7, 8 >v 


X (2) 3 ' 4 




\ <S) 1,2 






(7) 13, 14 ^ 


V (5) 9, 10 


Dipole 4 — *■ — S» 


(6) 11, 12 s 


Dipole 2 



Dipole 3 



Figure 3: The four categories of dipoles (Dipole 1, 4), the seven possible splittings and their order of 
creation in AutoDipole, see Table [2] for details. 



1. Choose all possible emitter pairs from the external legs. 

For a given real emission 2 — » (n + l)-particle scattering process we abbreviate the set of initial 
(final) state partons collectively as {initial} ({final}), i.e. the scattering reaction reads genetically 

{initial} -» {final} . (3.1) 

Then, the first step in the construction of the dipoles is the choice of the emitter pair, that is the 
root of the splitting of the quarks and gluons. This also specifies the kind of splitting. At NLO 
in QCD, there exist seven possible splittings, which we group into four classes and we enumerate 
these types of dipoles accordingly, Dipole 1, and Dipole 4. The ordering is shown in Figure[3] 
Due to the factorization of the splitting each dipole is associated to a certain reduced Born process 
and the following equations show schematically the operations performed on the sets of partons in 



Eq. (|3.1I) in order to construct the corresponding reduced Born process, 

Dipole 1: {initial} -> {final} -g, (3.2) 

Dipole 2: {initial} -» {final} -ff + g, (3.3) 

Dipole3: {initial} -/ + g {final}-/, (3.4) 

Dipole4: {initial} -g + f -> {final}-/, (3.5) 



where g is a gluon and /(/) a quark (anti-quark). The notation in Eqs. (13 .21) — (13 .51) indicates which 
partons are removed from or added to the respective sets {initial} and {final}. The cases Dipole 3 
and 4 can also occur with a reduced Born cross section where / and / are exchanged. The groups 
Dipole 2, 3 and 4 also need to account for the existence of various quark flavors, i.e. they exhibit 
manifest flavor dependence, see Table |2] for details. 
2. Choose all possible spectators for each emitter pair. 

The spectator is one external field which is different from both fields of the emitter pair and its 
choice is a purely combinatorial procedure. We use indices, i,j, and k, for fields in a final state 
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and, respectively, indices, a and b, for an initial state field. For a spectator in the final (initial) state 
denoted by k (b), this condition means k ± i,j (b + a). It emerges from a special feature of the 
subtraction formalism namely that the squared matrix elements with the Casimir operator, ( ), 
is expressed, due to color conservation, through CLBS < T, • )(i + k) which always includes one 
color dipole configuration in the part. 

3. Construct the dipole terms from the chosen combinations of emitter and spectator. 
The previous steps provide all such combinations as pairings of the type (emitter, spectator)= 
(ij,k),(ij,b),(ai,k), and (ai,b). Each case corresponds to one dipole term, for which we use the 
short-hands D^, D?., D^', and D ai ' , respectively. Detailed information about all pairings of 
emitter and spectator as well as the potential presence of a mass parameter is given in Table [2] 
in Appendix lAl For explicit expressions we refer to [EDS- 

Here it is worth to stress a few points. First of all, the chosen order for the generation has the 
advantage of grouping together the reduced Born matrix elements. The CLBS in the same category 
exhibit the same Lorentz structure, only the extended color matrix is different. This ordering leads 
to a block-diagonal structure in the evaluation of the color correlations. It leads to flat Fortran code 
including higher readability. Additionally, for the category Dipole 1 the Casimir operator Ty. in the 
denominator always cancels against the same one in the dipole splitting function. This cancellation 
leads to manifest simplifications. Finally, in our set-up we have access to the symbolic expression 
and can easily identify or extract partial subsets of dipole terms, if needed. For example we may 
want to discard the r-F(heavy-quark) splitting in the category Dipole 2, because this splitting does 
not give rise to poles in e. Our way of generating the dipoles in the package easily allows this 
selection. The user can discard dipoles of a specific category in the input file parameter . in (see 
Section HI). 

At this stage it remains to discuss the generation of all terms involving the I-, P- and K-operators. 
This is sketched in Figure @J Given, that these operators originate from the phase space integral 
over unresolved parton of the dipole terms, it is clear that we only need a subset of all previously 
generated terms. From all reduced Born amplitudes generated for the regularization of the real 
emission contributions (i.e. Figure [3]) we need only the CLBS in the category Dipole 1 for all 
I-, P-, and K-operator insertions. Hence, no new information is required and the Mathematica 
package assembles the respective expressions (along with the color correlations, see Section [331 ) 
in the Fortran files I term. £ and PKterm. £. 

3.3 The color linked Born amplitude squared 

The numerical evaluation of Born amplitudes (or their squares) is a routine task for many publicly 
available packages designed for automated LO calculations. As announced above, we choose 
the stand-alone version of MadGraph |f3]-[5l for this purpose and the added feature needed is the 
evaluation of the color link operators, see e.g. Eq. (12.71) . This can be done by a patch and allows 
us to obtain all CLBS in an automatic way, see Figured 

The generation of the amplitudes in MadGraph proceeds via Feynman diagrams and the color 
factors are separated from each diagram. During the evaluation everything is expressed in terms 
of generators of the fundamental representation of the color SU(3). A typical example is that the 
factor f abc of the gluon three-point vertex is rewritten in terms of the fundamental generator t a . 
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(6) 11, 12 

Figure 4: The order of creation of the integrated dipoles in AutoDipole with insertion of the I-operator 
(left) and of the P- and K-operators (right), see Tables[3]and |4]for details. 



with the help of the identity, 

f abc = _ 2 - ( Tr[ f t b f] _ ^ x \t c t b t a ] ) . (3.6) 

The color factors of each diagram are sorted in a unique order and they are expressed in a sum. 
When a specific term of a diagram is identical to one of the other diagrams, both are combined as 

M = ^C fl J fl , (3.7) 

a 

where C a denotes the independent color factors. Each C a has fundamental and adjoint color indices 
corresponding to the external quarks and the gluons, respectively. J a is the joint amplitude, e.g. 
Ji = +A\ - A3 H — where A,- is the partial amplitude of i-th diagram (with the color factor stripped 
off). The invariant matrix element squared is finally expressed in the form, 

|M| 2 = CF J, (3.8) 

where the color matrix CF is defined as 

(CFU = J] C*C b . (3.9) 

color 

For the CLBS we need to evaluate Eq. (13.81 ) with an insertion of two additional color operators to 
the emitter and spectator legs. This is precisely what our patch of MadGraph does (see Figured] 
and the directory lib/patch). The subroutines of MadGraph for the color factor calculations are 
well structured and the original routines to add the color factors t?. and f abc can be applied to the 
additional color insertions for the CLBS. The color algebra of the SU(3) is performed numerically 
and the resulting extended color matrix CF is written to the file allcolormat . inc, see Figure[2] 
We have realized the two color operator insertions for all CLBS in an automatic way and we have 
also checked that MadGraph with our interface works for a set of rather involved processes. As a 
simple example let us discuss the color insertions required for the process g(a)g(b) — > w(l)w(2)g(3). 
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The reduced Born process g(a)g(b) — » w(l)w(2) has three diagrams and the color factors are com- 
bined into two independent ones, (Ci,C 2 ) = {(t a t b )\2,{t b t a )n)- The components of the color matrix 
are given by the traces, (CF) n = (CF) 22 = Tr[t b t a t a t b ] and (CF) 12 = (CF) 2 i = Tr[t b t a t b t a ]. Then the 
color matrix is calculated as 

CF = ('^:^l (3.10) 



-2/3 16/3 

In the CLBS we need for instance the fundamental operator insertions into the legs 1 and 2. 
The components of the color matrix CF are modified to (CF')n = Tr[t b t a ft a t b t c ] and (CF')i2 = 
Tx\t b t a ft b t a f]. Then the modified color matrix is obtained as 

One of the more complicated examples consists of the two color operator insertions into the process 
g(a)g(b) — » f(l)?(2)g(3)g(4). In MadGraph the normal SU(3) color matrix for the process is a 
24 x 24 matrix. The first 15 components in the first row read (we refrain from spelling out the rest), 

CF = ^-(512,8,-64,80,8,-10,-1, -64, -64,8,-1, -10,-1, 62,-10,---). (3.12) 

Upon insertions of two adjoint color operators for a gluon into the legs 3 and 4 the extended 
routines calculate the modified color matrix as 

CF' = 1(8,0,8,16,0, -2,0,8,-1, -1,1, 2,-8,-7,1, •••), (3.13) 

and, of course, the result in Eq. (13.131) agrees with independent checks. 



3.4 Spin correlations for gluon emitters 

As we have seen above the dipoles for the splittings g — » gg and g — » qq (which involve a gluon 
emitter) introduce spin correlations (see e.g. Eq. (12.71) ). Since our numerical evaluation of the 
reduced Born amplitudes uses helicity amplitudes, we have to derive the components of the tensor 
for the spin correlations in the helicity formalism as well. To that end, we choose the definitions 
in [|37ll which we here call the XZC gluon polarization vector. 

Let us illustrate the necessary steps (formulated in 113810 with the case of the massless dipole term 
D(gi, gj ,k) inEq. (O, 

D( 8i ,gj,k) = -—L.^ij-jAjTt.TyAv), (3.14) 



where V^ v g k is written schematically as 



V^ g . k = l6na s M 2£ C A (-C lg ^ + C2^L v ), (3.15) 



in 



with If = HP; -ZjPj and some Lorentz scalars Ci ;2 (e.g. given in Eq. (12.81) ). Moreover. 
Eq. (13.141) we have expressed the CLBS in Eq. (12.91) through amplitudes and A v , where the 



12 



polarization vector e^(pij) of the emitter gluon with momentum pij has been amputated. Then we 
can transform both the dipole splitting function V^ v and the amplitude A„ to a helicity basis by 
inserting the polarization sum Y,a e fi* e v as 

A;V» V A V = YjA* a ,V a ' a A a , (3.16) 

A', A 

v a'a = ^'*y^y^ t (3 . 17) 
and we obtain the dipole D(gj,gj,k) in a +-helicity basis as 

D(gi,gj,k) = -^—I6na s ^ (3.18) 
3-Pi-Pj 

x [(d + C 2 \E + \ 2 ) (A* + T k ■ T,jA + +A*_T k • TyA_) + 2 C 2 Rc(E* + E.A* + T k ■ T ;J A_)] . 

Here we have used gauge invariance, the on-shell condition for the gluon and orthogonality of 
the polarization vectors, which gives rise to the relations e A * ■ pij = A ■ pij = L- pij = p 2 . = 0. The 
quantity E+ needed to express the dipoles in a helicity basis is defined as (see e.g. 113811 ) 

E ± = 6+-L. (3.19) 

As a subtlety, we would like to point out the following. In the HELAS library (6l|Vl the gluon 
polarization vector is calculated by the subroutine VXXXXX and taken to be in a helicity basis but 
the phase conventions are different from the ones in OTTl . Thus we have to relate these conventions 
to our choice |[3~7Tl . Although it may be natural to use the same definition for the gluon polariza- 



tion vectors as in HELAS for the dipole splitting functions (e.g. V^ g . k m Eq. (13.151 )) we have 
the freedom not to do so because of gauge invariance and the on-shell condition. In our imple- 
mentation we have chosen the latter option with the advantage that the calculation of the dipole 
splitting functions can be completely separated from the part for evaluation of the CLBS. With the 
XZC definitions E+ in Eq. (13.191 ) is expressible in terms of spinor products which can be easily 
implemented in Fortran code. 

The polarization vectors of HELAS [6,7] are written in a helicity basis as 

<[HELAS](M) = -^( + ^(D- Z ' e M2))' (3-20) 

with the explicit expression for the vectors in a linear basis, 6^,(1,2) given in (6113. Likewise, for 
XZC W\ we have 



{ q-\jfi\k- ) 

^J2{qk)* 



where the relation e + r * 1 = e r , holds and q is an arbitrary reference momentum. The difference 

^/[XZCJ fi\. XZC I 

between the conventions Eq. (13.201 ) and Eq. (13.211) amounts to a complex phase, 



«£]<toHW*'90 = (3-22) 

The phase difference contributes to only the second term in the square bracket in Eq. (13.181 ) and it 
is then rewritten as 

+ 2C 2 Re(E;E_A;T,-T ;j A_) [xzc] = -2 C 2 Re((e + '^ E+ f(A* + T k • TyA-)^) , (3.23) 
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where the CLBS is also taken in the helicity basis for the gluon emitter as in HELAS. The quantity 
E+ = j • L in Eq. (13.231) can be computed according to Eq. (13.191 ) in the conventions of 113711 in 
terms of spinor products as (see e.g. [1381 ) 

e +W.k) E+ = e +WPij) z J^EJPA tjyPjf (3 24) 

^(Pjpij) 

where we have set k = pij and chosen q = pj for the reference momentum in Eq. (13 .2 1 b . Explicit 
expressions for the vector E+ in Eq. (13.191 ) for all required momentum configurations (massless 
and massive) are given in Appendix [Bj For the dipole terms including the massive partons, the 
result for E+ in Eq. (13.191 ) contains terms like (p, - \pi\Pai ~ ) with on-shell momenta p 2 = m 2 . In 
order to express these terms through spinor products we have used flat momenta, 

2 

171 

p\ = Pk-~ n k —Pai, (3-25) 

2pk-Pail 

which are massless because (p\) 2 = 0. Due to the equations of motion p a i \pm - ) = 0, we can 
rewrite the quantities under consideration, e.g. 

(Pi-\h\Pai-) = {Pi-\p\\Pai-), (3-26) 

so that they are accessible to standard spinor calculus, see also [1391 for the spinor helicity for- 
malism including massive fermions. Our treatment of the spin correlations is contained in the 
Mathematica sources where the expressions for E+ (as given in Appendix |B]) are implemented. 
During the automated generation of all subtraction terms for numerical evaluation the results are 
written to the output in the Fortran file, dipole . f , see Fig. [2] 



3.5 Checking the soft and collinear limits 

Let us discuss the checks for the AutoDipole package. First, there is the standard quality check 
on the automatic generation of the dipole terms because the result according to Eq. (12.31 ) has to be 
finite when approaching the singular regions. Second, since the Catani-Seymour formalism is well 
established, we can compare the output of AutoDipole with results of independent implementations 
in the literature. 

The AutoDipole package performs automatically for all generated subtraction terms a numerical 
check of all IR limits. For that purpose the subtracted squared matrix element for the process 
under consideration are constructed and sampled over all limits to test the cancellation of the IR 
singularity. For a given soft/collinear limit L ; , we pick up the set S(i) of the corresponding dipoles 
to test whether the quantity 



lim 



|M| 2 - £ D(j) 

jcS(0 



(3.27) 



is soft/collinear safe. The code lists all limits L ; , and the corresponding sets S(z). A concrete 
example will be presented in Subsection 14.21 The leading soft/collinear singularity of the squared 
matrix element |M| 2 behaves in the soft limit as 1 fk 2 for a gluon of momentum k and in the collinear 



14 



one as 1 /(2pi ■ pj) = 1 / Sij for two collinear momenta pi and pj. The set of the corresponding dipoles 
cancels this leading singularity pointwise as 

|M| 2 - ^ D(j) = —(a + aix + a 2 x 2 + ■■■)- ^a , (3.28) 

jcS(0 X X 

where x = k( -Jsjj) is the respective Lorentz scalar depending on the momenta in the soft (collinear) 
limit. The difference in Eq. (13.281) is integrable over the real emission phase space. The numer- 
ical accuracy, however, has to be controlled, because deep in the singular regions the individual 
contributions in Eq. (13.281) become very large and with the limited numerical precision of floating 
point arithmetic the cancellation might be imperfect. Fortunately, the stability of the numerical 
cancellation can be tested automatically. Normalizing Eq. (13.281 ) one can check the slope of the 
quantity, 

|M| 2 -ZD 1, M 2i , nw 
— — = —(aix + a 2 x +■■■). (3.29) 

]M] 2 a 

Based on Eq. (13.291) AutoDipole tests the cancellation and returns values for the fiducial regions 
for the Lorentz invariants Sjj. 

We have extensively tested the soft and collinear finiteness of the generated subtraction terms for 
numerous scattering processes in e + e~ , ep and pp-collisions including massive quarks and weak 
gauge bosons. Among those are a large class of 2 — » 4 scattering processes as needed in the NLO 
QCD correction to 2 — » 3 reactions. We have also tested scattering processes with 2 — > 5 and 2 — > 6 
partons, like e.g. 

uu — » ddggg 
gg -> W + udgg 
gg -» ttggg 
gg -> ttbbg 
gg -> ttbbgg 

which are currently under investigation in view of phenomenological applications for the LHC. 
In addition, we have been able to obtain also perfect agreement with published results H36U401I . For 
the NLO QCD corrections to pp — » tt + ljet production results for the real emission contributions 
have been presented at individual phase space points in IT361 . We agree at least to 14 digits for the 
matrix elements squared and at least to 12 digits for the sum of the subtraction terms. All details 
of our comparison are shown in Table [5] in Appendix O We have also checked against the results 
for the NLO QCD corrections to pp — » + ljet [|40ll , in particular the real emissions in the 

channel uu — » W + W~gg. Also here we have obtained very good agreement as documented in the 
last entry of Table [51 The AutoDipole package provides a shell script to reproduce some of the 
numbers in that Table [5J (see Section 1441) . 

Let us briefly comment on the integrated dipoles in Eq. (12.31) originating from the I-operators. 
Here, a complete check is more involved because for a given process scheme dependence enters, 
i.e. whether the singularities of the NLO virtual contributions are factorized with respect to the 
J-dimensional or 4-dimensional Born amplitude. 
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Nevertheless, partial checks of the singularity structure are straightforward. The matrix element 
squared with the sum of all I-insertions is written schematically as 

<l,...,n|I|l,...,n> = C-2^ + C_i- + C , (3.30) 

where C-\ and C_2 are process (and kinematics) dependent coefficients. For the leading pole 
~ 1 /e 2 , the coefficient C_2 obeys the following simple relation 

C-2 = ^\M B J 2 ((n q + nij )C F +n g CA), (3.31) 

where M Bmn is the corresponding Born amplitude. The number of external massless (anti-)quarks 
is n q (n q ) and n g the number of external gluons. In QCD the standard SU(A0 color factors take the 
numerical values Cf = 4/3 and Ca = 3. 

Apart from simple relations like Eq. (13.311) the coefficients for color and spin averaged results for 
the I-operator could again be compared with the published literature for pp — » tt + ljet production 
at NLO in QCD [36J. We have obtained agreement to at least 14 digits as shown in Table [6] The 
numbers in Table [6] can again be reproduced with a simple command explained in Section I4~4l 



4 Running AutoDipole 

This Section is meant to be a short manual to the AutoDipole package. We explain the installation 
procedure, the use and discuss a few processes as examples. 

4.1 Install 

The package, AutoDipole_Vl .2.3. tar, is available for download from PTTl (or else from the 
authors upon the request). In addition, one must also obtain the stand-alone version of MadGraph 
(MG_ME_SA_V4.4.30) e.g. by download from 021 and put it (in the form of . tar . gz or . tar) 
in the AutoDipole package directory. Then execute the installation procedure by3 

./install . sh 

and the directory structure as displayed in Figure [Demerges. Repetition of this command always 
allows the user to recover the initial settings described in this Section. 

4.2 A short example 

Let us next illustrate the use of AutoDipole. First of all, one has to specify all partonic real emission 
processes which appear at NLO in the observable under consideration. To that end, as a concrete 
example let us choose a simple process, uu — > ddg which contributes to hadronic di-jet production, 

1 The AutoDipole package (version 1.2.3 and earlier) is assumed to be used with MadGraph (stand-alone version 
MG_ME_SA_V4.4.30). Upgrades to newer versions of MadGraph require the user to set the MadGraph version in the 
first line of install . sh by hand before installation. We have tested that AutoDipole version 1.2.3 also works with 
the MadGraph stand-alone version MG_ME_SA_V4.4.39. An automated test of the installation can be executed with 
the command check . sh as detailed in Section l4~4l 
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pp —> 2 jets at NLO. This example exhibits all features of the Catani-Seymour formalism, i.e. it 
needs dipoles and I-operators as well as the P- and K-terms. We start with the Mathematica part 
of AutoDipole. The package can simply be included through the driver file as 

«driver_user . m 

Next, we can run the package for the process uu — > ddg with the command 
GenerateAll[{u,ubar} , {d,dbar,g}] 

Upon running the AutoDipole package successfully, all generated Fortran files are stored in a 
(newly created) directory under process, for the example at hand . /process/Proc_uux_ddxg. 
AutoDipole returns the message: 

Run has been succeeded. 

In order to run the generated code, please enter 
in the newly created directory under ./process 
make 
. /runD 



The real matrix element squared |M| 2 and the sum of all dipole terms 2,-DO") are evaluated at 10 
phase space points with this command. Also the check of cancellations in all infrared limits for 
the subtracted matrix element can performed. The respective command is: 

make checkIR 
./checkIR 

Likewise, the integrated dipoles (I-operators) as well as the contribution of the P- and K-terms 
at 10 points in phase space can be evaluated in the subdirectories Proc_uux_ddxg/Virtual and 
Proc_uux_ddxg/PK with the commands: 

make runl 
. /runl 

make runPK 
. /runPK . 

The present example uii — > ddg runs with the default settings of AutoDipole. For further illustration 
the package also comes along with a prepared list of examples so that explicit numbers can be 
obtained by executing the shell scripts explained in the next subsection. 

If the user wants to include the generated Fortran code with the subtraction terms into his/her own 
project for further numerical evaluations, he/she needs (i.e. has to copy) a number of files for that 
purpose. The complete list of generated files is given in Table \T\ Note that most of the functions 
are only needed internally and are not meant to be called by the user directly. We assume here that 
all input momenta {/?,} are generated from phase space routines supplied by the user and that the 
HELAS library libdhelas3 . a is linked. 

Dipole terms : directory process/Proc_uux_ddxg 

From the file check . f the call of the relevant subroutines for the evaluation of the dipole terms is 
obvious. The user needs to copy the following files: 
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file (subroutine) 



input 



output 



[Dipole terms] 
matrix, f (smatrix) 
dipole . f (dipole) 
reducedm.f (reducedm) 
cmatrix. £ (cmatrix#) 
cmatrix. f (cmatrix#dh) 
check . f (program checkdipole) 

[I-terms] 

matrixLOl . £ (smatrix) 
I term, f (Iterm) 
cmatrix. £ (cmatrix#) 
checkVirt . £ (program check Virt) 

[P- and K-terms] 
matrixLOl . £ (smatrix) 
PKterm.f (PKterm) 
cmatrix. £ (cmatrix#) 
checkPK . £ (program checkPK) 

[Check of soft/collinear limits] 
collinear . £ 

checklR. £ (program checkIR) 



[Pi) 
iPi) 
{Pi\ 

reduced [pi\ 
reduced {pi} 
{pi} in inputm.h 



{Pi) 
{pi) 
{Pi\ 

{pi} in inputmLO.h 



{Pi) 
iPi) 
{Pi} 

{pi} in inputmLO.h 



{Pi) 

collinear limit of {p{\ 



|M| 2 for real emission 
D(0 

reduced {pi} 
CLBS (like helicity) 
CLBS (unlike helicity) 
check of D(i), |M| 2 



|M Bora | 2 

Ci of I-terms 

CLBS 

check of Ci of I-terms 



|M Bora | 2 

Pt(z), K(0 (P- and K-terms) 
CLBS 

check of Pt(i) and K(i) 



collinear limits of 
(|M| 2 -ZD)/|M| 2 



common input files 



[Process dependent] 
allcolormat# . inc 
emitspecinfo . inc 
param_card.dat 
nexternal . inc 

nexternal2 . inc 
inputm.h 
inputmLO . h 

[Process independent] 
coupl . inc 
couplings . £ 



extended color matrices (used in cmatrix#) 
information of emitter for dipoles (used in cmatrix#) 
parameters for CLBS 

field numbers of real emission process for dipoles and 

of Born process for I-, P- and K-terms 

field number of reduced Born process 

phase space points for dipole checks 

phase space points for I-, P-, and K-term checks 



common block of parameters 

reads parameters in param_card . dat and assigns 

them in coupl . inc 



Table 1 : The file content of the AutoDipole package. In names of subroutine and files (e.g. cmatrix#) the 
symbol # denotes an integer number, typically #=1. 
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matrix. £ 
dipole . f 
reducedm. f 
cmatrix . £ 
allcolormat\# . inc 
emitspecinfo . inc 
nexternal . inc 
nexternal2 . inc 

coupl . inc 
couplings. £ 

param_card.dat 



(generated by AutoDipole) 



(generated by AutoDipole) 

(process independent file) 
(process independent file) 

(parameter input file) 



To use the generated code, one first has to initialize MadGraph. The parameters are determined in 
param_card . dat and are assigned to the definitions in coupl . inc with the help of the subroutine 
setpara in couplings . f as 



include 'coupl. inc' 

call setparaCparam_card.dat' 



.true.) 



to be used subsequently in the subroutines cmatrix#, cmatrix#dh, and smatrix through the 
common blocks. The symbol # here and above denotes an integer number (e.g. #=1). For the 
dipoles only the strong coupling constant a s as well as the top and bottom masses are needed in 
the evaluation. These three parameters appear as 

AL=g**2*(4.dQ*pi)**(-D 

mt=tmass 

mb=bmass 

and are passed to the subroutine dipole through the common blocks, 

common /usedalpha/ AL 
common /MASS/ mt,mb 

Next, the evaluation of the squared matrix element for the real emission process |M| 2 proceeds in 
exactly the same way as in MadGraph with the call 

call smatrix(p,msq) 

Finally, the dipole terms are evaluated by calling the subroutine 
call dipole (p, dip, SumD) 

which returns an array with the values for each dipole term (dip) and the sum of the all dipoles 
(SumD) as output. For our example process uu —* ddg we have 15 individual dipoles. The reduced 
kinematics of each dipole (as determined by the subroutine reducedm) can be accessed by the 
common block 
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double precision ptil(0: 3 , 1 :nl ,n2) 
common /OUTPUTM/ ptil 

where nl=4 and n2 = 15 in the definition of the momenta for the example at hand. 
I-terms : directory process/Proc_uux_ddxg/Virtual 

Our program checkVirt . f for checks of the I-terms again provides an example for the use 
of the subroutines to evaluate the I-terms. In addition to the eight files from cmatrix.f to 
param_card . dat listed above for the evaluation of the dipoles, the user has to copy one file: 

Iterm. f 

Then, the initialization phase is exactly the same as before and for the evaluation of the I-terms, 
one has to call the subroutine 

call Iterm(p,coe£,SumI) 

which provides the value of the coefficients (i.e. the array coef) as well as the sum (SumI) of all 
I-terms as output. If the user also wants to compute the LO matrix element |M Bom | 2 he needs the 
file matrixLOl . £. The execution of the subroutine with the command 

call smatrix(p,msqLO) 

is performed again exactly as in MadGraph. 

P- and K-terms : directory process/Proc_uux_ddxg/PK 
The user needs the following file to evaluate the P- and K-terms: 

PKterm. f 
matrixLOl . f 

as well as the eight files from cmatrix . £ to param_card . dat already discussed for the dipoles. 
The initialization phase is unchanged but the user has to provide the parton momentum fraction x 
from the mass factorization as additional input. The P- and K-terms for the first leg are computed 
by the call 

call PKterml(p,x,SumP,SumK) 

which returns separately the sum of all P- and K-terms (SumP and SumK). For the second parton 
in initial state (as in hadron-hadron collisions) the corresponding P- and K-terms are evaluated by 
calling PKterm2. 

4.3 General usage 

For a specific application of AutoDipole, the user needs to perform the following three steps: 

1. Setup of parameters : parameter. m 

2. Run of package : GenerateAll [{initial} , {final}] 

3. Run and checks of the generated code at : /process/Proc_xx_xx/ . 
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Let us explain each of these steps in more detail. 
1. Setup of parameters 

Generally, the user has to supply all process parameters needed for the dipole subtraction procedure 
by editing the file, 

parameter .m 

where (following the MadGraph conventions) all values with the dimension of a mass are in units 
of [GeV]. Let us explain all variables in the file together with the chosen default values. 

Parameters for dipole terms 

ep=0 : 

The parameter of dimensional regularization of space-time, D = 4 - 2e. The default choice 
eliminates higher orders (positive powers) of e in the dipole terms. 

kap=2/3 : 

k is free parameter in some expressions for massive dipoles 0. The value k = leads to the 
simplest results for the dipole terms, while k = 2/3 produces the simplest expressions of the 
I-terms. 

skipdipole={ } : 

This set specifies the kind of dipoles to be skipped during the creation. For example, the set 
skipdipole={2t} omits the creation of the Dipole 2-(5) with the t — t splitting (see Table [2] 
in Appendix lAl). 

mur= 174.30 : 

The value of the renormalization scale. The default value is the same as the default value of 
the top-quark mass in MadGraph. 

acccut=10 A (-3) : 

This value is used for the checking the quality of the infrared cancellations in the generated 
code upon scanning over all possible collinear configurations. 

Parameters for I-terms 

mFlist={ } : 

This set lists all heavy quarks which possibly contribute to the I-terms with a gluon emitter 
in final state, e.g. as mFlist={t,b}. 

lightflavors=l : 

This is the number of the light flavors which contributes to the I-terms with a gluon emitter 
in initial and final states. 

replistVirtual={ } : 

This set enables user defined symbolic replacements to simplify the virtual contribution, e.g. 
of the type, replistVirtual={Gamma[l-eps]-> (4*Pi) A (eps)*somesymbol A (-l)} which 
results in the replacement, Xrf . =somesymbol. 
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si jexpansion=0 : 

This parameter determines whether the typical factors of dimensional regularization like 
j in the I-terms are expanded (0) or not (1). 

Parameters for P- and K-terms 

KFSf f =0, KFSgg=0, KFSfg=0, and KFSgf=0 : 

Theses parameters specify the factorization scheme. In the MS scheme (default choice) they 
are all vanishing. 

lightquarksPK3={u} : 

This set specifies all light flavors to be added to the initial gluon of the reduced Born process 
in the P- and K-terms (see group (6) in Table |4] in Appendix lAl). 

muFfort= 174.30 : 

This is the value of the factorization scale /Uf , which has the same default value as the 
renormalization scale /j.r. 

The other parameters which are used by MadGraph are determined in the file 
patch/Models/sm/param_card . dat. 

The default values for the most important ones are the following (in units of [GeV] for masses), 
a,=0.1180 

aw=0.00754716981 1320755 

Top mass=174.3 

Bottom mass=4.7 

W mass=80.419 

Zmass=91.188 

Higgs mass= 120.0 

These parameters are used consistently throughout the run of the package. Any other changes of 
SM parameters like e.g. Gf, the CKM -parameters, the decay widths of heavy fields, and so on, are 
done in the same MadGraph file, param_card . dat. Similarly, also control over all interactions in 
the SM proceeds entirely through MadGraph, i.e. through the file, 
patch/Models/sm/interactions . dat 

in the exactly same way as in the normal use of MadGraph. All MadGraph parameters are ex- 
plained in ma. 

2. Run of package 

For general QCD processes (including leptons and SM bosons) we run AutoDipole by executing 
the Mathematica command: 

GenerateAll [{initial} , {final}] 

where the sets labeled initial and final contain all external particles of the respective scattering 
process. The following fields are available in the present version of the package: 
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Fermions: 

(u,d,b,t,e~ ,ju~ ,u,d,b,t,e + = 

(u, d, b, t, e, muon, ubar, dbar, bbar, tbar, ebar, muonbar) , 

Bosons: 

(gluon, photon, W + , W~, Z, Higg s) = (g, gamma, Wp, Wm, Z, h) , 

where the right hand side denotes the input format of AutoDipole. The available quark fields 
include the first family (u and J-quarks) which serves as a template for light quarks and, respec- 
tively, the third family (b and ^-quarks) as one for heavy quarks. In practice this suffices as far as 
the subtraction formalism is concerned, because generally QCD scattering processes at NLO will 
be insensitive to the individual flavor of the quarks, and the sum over all light quark flavors can be 
trivially taken afterwards. If one intends to study particular channels, e.g. in W ± production, flavor 
separation is required, though. 

As a result of the run, AutoDipole produces the following information about the created dipoles 
and the I-, P- and K-terms: 

- All analytical expressions for the dipole terms are written in the order given in Table [2] along 
with the reduced Born process (emitter, spectator) and the reduced kinematics for each dipole 
term. Also, the number of created dipole terms is shown group by group, namely Dipole 1, 
2, 3, and 4 and their subgroups (see Figure©. The number of massive dipole terms is given 
separately. 

- All soft and collinear limits together with the corresponding dipole terms are returned. 

For example, for the process g(l)g(2) — » ?(3)?(4)g(5)g(6) (see next Subsection) the fifth 
collinear limit is denoted Lc(5) = (5,6) : Sc(5) = {17, 18, 19,20}. This means that the collinear 
divergence of the pair (5(g), 6(g)) is canceled by the set of the dipoles, Sc(5) and the set 
consists of the dipoles number 17,18,19, and 20. In the same way, the soft limit and the cor- 
responding dipoles are displayed. E.g. for the soft limit of gluon 5(g), Ls(l) = (5) : Ss(l) = 
{Sc(l),Sc(2),Sc(5),Sc(6),Sc(7)}, which means that the soft divergence of 5(g) is canceled 
by the set of the dipoles, Ss(l). This set in turn consists of the sets of collinear limits Sc(l), 
Sc(2), and so on. 

- The analytical expressions of all I-terms are shown in the order as given in Table |3] prior to 
the expansion in e as well as the individual coefficients of the poles C_2, C_i, and Co (see 
Eq. (|3.30l) ) after the e expansion. Also the total number of all created I-terms is displayed 
by group by group (see Table [3]). 

- The analytical expressions of all created P- and K-terms are shown in the order given in 
Table |4] together with the total number of P- and K-terms. 

- Finally, AutoDipole outputs a summary with the numbers of all created dipoles, I-, P- and 
K-terms and the reduced Born processes. 

This information, being rather detailed, helps with the identification of individual terms according 
to the order of creation. The Mathematica part of AutoDipole provides analytical and symbolic 
expressions, which are easily accessible, if needed. These features are primarily for an experienced 
user. 



23 



Upon running the AutoDipole package it occurs (on purely algorithmic grounds) in some cases 
though, that a particular reduced Born cross section does not exist. For example, the process 
e~e + — » uug does give rise through the dipole 2-(5)-(w,w) to the Born process e~e + — » gg, which 
does not exist at tree level. In such cases, an error message is displayed: 

Reduced Born amplitude B2u does not exist. 
Dipole 2u must be switched off. 
Please set in the file . /parameter. m 
skipdipole={ 2u } 

Following the message the user has to modify the command skipdipole and repeat the Math- 
ematica run. In this way only the relevant partial subsets of dipole terms is extracted and the 
offending Born processes are skipped (thanks to skipdipole) before calling the MadGraph part. 

3. Run and checks of the generated code 

Finally, as briefly explained above, the user can run the generated Fortran code for the regularized 
real emission process as well as the I-, P- and K-terms. The checks proceed as follows: 

Dipole terms 

As mentioned above, the evaluation of the real matrix element squared |M| 2 and the sum 
of all dipole terms ZjD(i) at 10 randomly chosen phase space points as well as the check 
of cancellations in all infrared limits for the subtracted matrix element are automatically 
executed: 



make 
. /runD 



make checkIR 
./checkIR 



The latter command performs a scan over all possible collinear configurations. For the ex- 
pression 

|M| 2 -yD 

■ — =— (4 1) 

|M| 2 ' { } 

cf. Eq. (13.291 ), it is checked whether the cancellation works at least to the accuracy defined 
by the value of the parameter acccut in the file parameter. m. If the absolute value of 
the normalized subtracted squared matrix element in Eq. (13.291 ) drops below the pre-defined 
value of acccut sufficiently deep in the collinear limit, the message, 

Infrared safety of all collinear limits is confirmed for S_ij/S > ... 

is printed. 
I-terms 

The Fortran code to evaluate all I-terms is collected in the subdirectory ./Virtual. There, 
the coefficients of the double and single poles in the e expansion C_2> C-b Q) m Eq. (13.301 ) 
can be evaluated at 10 randomly chosen phase space points by: 
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make runl 
. /runl 

This includes in particular the calculation of the double poles according to 



(a s /2n)\M Bom \ 2 
and the user can easily confirm that Eq. (|3.31l) holds. 
P- and K-terms 

Likewise, the mass factorization terms originating from the P and K-terms is stored in . /PK 
and also evaluated at 10 randomly chosen phase space points by the command: 

make runPK 
. /runPK 

This concludes our brief description of the user specific features of the AutoDipole package. 
4.4 Examples 

The package includes a command to demonstrate the run for several example processes. This helps 
to explain the package and at the same time confirms that the installation was done successfully 
and that the package works properly. These examples include the computation of some of the 
numbers given in Tables [5] and [6] in Appendix O The command 

. /lib/sh/check . sh 

generates five processes in different kinematics, 

1. e~e + — » uug 

2. e~u — » e~ug 
3- gg^ttgg 

4. ug — » ttug 

5. uu^ W + W~gg 

provides suitable settings for the parameters of each process and runs the package. For illustration, 
we only explain in the following those settings which are changed with respect to the default ones. 
The processes from 1 to 4 use a value of a s =0. 1075205492734706 from ll36l for the strong cou- 
pling constant. For the processes 1 and 2 the contribution of the Z-boson exchange is switched off 
in the file interactions . dat for simplicity. Moreover, in parameter .m the splitting 2u in pro- 
cess 1 (and 3u in 2) is skipped by the command skipdipole = {2u} (and skipdipole ={3u} 
in 2, respectively) to exclude an inexistent reduced Born process. 

For the processes 3 and 4 the QED interaction is switched off in the file interactions.dat 
in order to focus on the QCD corrections only. Furthermore, a value of 174 [GeV] for the top 
mass is used following ll36l and the decay width of the top-quark must be switched off in the file 
param_card.dat. For the comparison with 11361 the t-t splitting is skipped in parameter. m, 
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i.e. skipdipole ={2t}, and the renormalization and factorization scales are set equal to the top 
mass value. For the I-terms five light flavors are taken into account, lightflavors = 5. 
Finally, process 5 uses the following settings PfOll in the file param_card.dat: W-mass=80.425, 
a s = 0.1202629003906369, a em =0.007543595708669335 and Fermi's constant G F = 1.16637 x 
10~ 5 , and the decay width of W-boson is switched off. The results of processes 3, 4, and 5 should 
reproduce the corresponding numbers in Tables |5] and |6] 

The test run of all examples is finally checked by comparing all computed values of the dipoles, 
I-terms, and the averaged matrix element squared |M| 2 for the processes 1 to 5 the with ones in 
the data files in . /lib/check/. At the end of the run, e.g. for example 3, the following message 
appears: 

3.gg_ttxgg 

|M| A 2(Real) : Confirmed 

Dipole : Confirmed 

|M| A 2(L0) : Confirmed 

I- term : Confirmed 

This confirms the successful installation and proper functioning of the package. The generated 
code for all examples is stored in the directory . /process as in a usual run and the directory of 
each process includes in the file res_std the result of the comparison with the data files from 
. /lib/check/. All parameters used for each process are accessible in param_card . dat, so the 
user can also run the respective files individually with usual commands of the previous subsection. 



5 Conclusion 

We have presented the package AutoDipole which provides an implementation of the Catani- 
Seymour dipole formalism (US to compute subtracted matrix elements squared in an automatic 
way. The package (partial aspects of which have been discussed before |fT71[T8l ) consists of es- 
sentially three ingredients: the automatic generation of all dipole terms via Mathematica routines 
as well as the calculation of the color-correlated matrix elements and the evaluation of different 
helicity amplitudes with the help of MadGraph (SHU (i.e. a patched stand-alone version). 
We have presented a complete discussion of the set-up of the program and the details of our imple- 
mentation as well as the numerical evaluation of the squared matrix elements. Particular emphasis 
has been put on explaining the use of AutoDipole which has been illustrated with sufficiently non- 
trivial examples. The finiteness in the soft and collinear limits has been demonstrated for a large 
number of processes in e + e~, ep and pp-scattering. Also, checks for the complete set of dipoles 
have been performed by comparing to various processes in the recent literature. 
In comparison with other available software to generate subtraction terms for the real emission 
contributions AutoDipole has a number of advantages. First of all, the algorithm is implemented 
in a way so that analytic expressions are generated which are accessible to symbolic manipulation. 
Then, the chosen order for the dipole generation produces flat Fortran code and allows to handle 
scattering processes with 2 — > 5 and 2 — > 6 partons of relevance for LHC phenomenology including 
the computation of the CLBS. Finally, the structure of the package is highly modular and flexible 
for further extensions. 
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Future improvements of the AutoDipole package will address restricting the dipoles to small re- 
gions of phase space, as realized e.g. through cr-parameters 11431 . This is advantageous when 
performing the Monte Carlo integration over the phase space of the real emission process. The 
implementation of such cuts affects also the integrated dipoles and the case of massless partons 
is available in the literature whereas for the case of massive partons most cases can be found 
in Il44ll45l . The remaining integrals for exhaustive coverage can be done e.g. with the methods 
of ||32l . Next, the current version of AutoDipole does not cover processes with identified hadrons 
in the final state where one needs fragmentation. Here, the so-called H-terms have to be taken into 
account, which will be implemented in a later version of AutoDipole. Finally, the structure of our 
code is particularly suited for extensions to physics beyond the Standard Model (BSM). The main 
modification needed is the so-called model file of MadGraph, which is currently already available 
for a large class of BSM scenarios including e.g. the minimal supersymmetric extension of Stan- 
dard Model for which the necessary modifications of the dipole formalism are known [2J. All the 
extensions mentioned above can be easily realized within the existing structure of our code. 

The AutoDipole package is available for download from PTTl or from the authors upon request. 
The MadGraph package (stand-alone version) which includes the HELAS library may be obtained 
from 11421 
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Appendix 

A Creation order of subtraction terms 
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Table 2: The emitter pair is denoted (/, j) or (a,i) and, respectively, the spectator k(± i,j) or ^(^ a). D is a 
massless dipole and DM is a massive dipole. The equation numbers refer to HI for massless dipoles (D) and 
to for massive ones (DM). Dipole 2 has the other flavors, massless down and massive bottom quarks. 
Dipole 3 and 4 have also the other light flavors, / = u,d,d. 
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Table 3: The emitter is denoted i or a and, respectively, the spectator k{+ i) or b{£ a). I is the massless 
insertion operator and IM is a massive one. The equation numbers refer to HI for massless insertions (I) 
and to [2] for massive ones (IM). The set of the heavy-quark masses the running in the quark loop are 
abbreviated {m F }. 
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Table 4: The emitter is denoted a' and different from spectator b{+ a'). The massless case is abbreviated K 
and massive one KM. The equation numbers refer to HI for massless insertions (K) and to for massive 
ones (KM). For the cases Dipole 3 and 4, / can take the light flavors f -u, u,d,d and so on. 



31 



B Correlation functions of gluon emitter 



We use the following definitions of spinor products 0711 : 

<*ifc 2 > = ^qiqe*^ - ^kpqe +i ^ k2) ©1) 
= \£\±k z dB2) 



e 



+ic/>(k) _ 



\k±\ 

k ± = k x + ik y (E4) 

If the reduced momentum of the initial gluon emitter lies in z-axis with I = £/( 1,0,0, + 1) the spinor 
product are defined as (pip a i) z to avoid an ambiguity in the intermediate steps of the calculation, 

j -2jEiE 2 (Ei = h z and E 2 = ~h z ) =,„ 
^ /l/2)z - 1 2V£i^ (Bi = and £ 2 = l 2z ) m) 



The explicit forms of e +l ^ E+ in Eq. (13.231) are shown for all massless and massive dipoles with 
a gluon emitter. Here we write the quantity including the phase difference as = e +l ^ E+. 
The reduced momenta are ones used in the corresponding equations in (HE), which are shown in 
Table [2l The pij\ and p a i are massless and the pM^ are massive ones. 

Massless momenta 

[Dipole 1-(2>3(to* = 0) and Dipole 2-(5)-9u(m fe = 0)] 



[Dipole l-(2)-4 and Dipole 2-(5)-10u(m jt = 0)] 
[Dipole l-(4)-7(m fe = 0) and Dipole 3-(6)-llu(m^ = 0)] 
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[Dipole l-(4)-8 and Dipole 3-(6)-12u] 

E4' + = -""' r _»-»'^>'«"»>: (B1()) 
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Massive momenta 

[Dipole l-(2)-3((ra& * 0)) and Dipole 2-(5)-9u((m* # 0))] 



M; a ,^y^^ ®ii) 



[Dipole 2-(5)-9t] 



^2{p)p~Mij) 



with flat momentum 

[Dipole 2-(5)-10t] 



^2(p)p~Mij) 



with 



7 2pifj)-pMi 



[Dipole l-(4)-7((m* * 0)) and Dipole 3-(6)-llu((m fc * 0))] 



with 
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©12) 



/7t 

"" pM f/ (El 3) 



EM2' + = e +l<fl(pM J 7 ' ^ © 1 4) 



/71 

b "" pM f/ ©15) 



_ . (PiP k )(P k Pai) z 

EM3' + = k — — ©16) 

V2(l -Ui){pip a i) z 



2 

p b k = Pk-^ J ^Pai ©17) 
2pk'Pai 



C Comparison with the tt + ljet 





MGeV" 4 ] 


^otGeV- 4 ] 


g(Pa)g(Pb) -> t(pt)i(pf)g(Pc)8(Pd) 


AutoDipole 
Ref. HSU 


7.82039670869613 • 10" 1U 
7.82039670869604(1) • 10" 10 


1.02594003852407 • 10~ y 
1.02594003852396(2) -lO" 9 


q(p a )q(pb) -> t(pt)t(pt)g(p c )g(Pd) 


AutoDipole 
Ref. [[21 


1.12077211361620- 10 _1U 
1.12077211361619(0) • 10" 10 


1.22619016939900- 1(T 10 
1.22619016939909(1) • 10" 10 


q(Pa)g{Pb) -> t(p t )t(pt)g(pc)q(pd) 


AutoDipole 
Ref. [[21 


2.75641273146785 -lO" 11 
2.75641273146783(0)- 10" 11 


4.79768338384667 -lO -11 
4.79768338384625(3)- 10" 11 


q(Pa)g(Pb) -* t(p t )t(pi)q(Pc)g(Pd) 


AutoDipole 
Ref. 01 


3.46150168295956- 10- 11 
3.46150168295954(1)- 10" 11 


8.34555795894942- 10- 11 
8.34555795894963(2)- 10" 11 


g(Pa)g(Pb) -* t(pt)t(pt)q(p c )q(pd) 


AutoDipole 
Ref. [El 


1.21420520114780- lO -11 
1.21420520114779(0)- 10" 11 


2.13553289076589- lO" 11 
2.13553289076550(3)- 10" 11 


q(Pa)q(pb) -> t(pt)t(pi)q(p c )q(pd) 


AutoDipole 
Ref. (21 


5.13710959990068 • 10" 12 
5.13710959990064(1)- 10~ 12 


9.06330902408356 • 10" 12 
9.06330902408275(3)- 10~ 12 


U(p a )u(p b ) W + (p w+ )W (p w -)g(Pc)g(.Pd) 


AutoDipole 
Ref. frfoTl 


0.627402537098012 -10~ 9 
0.627402537098007 • 10~ 9 


0.114149934878320- 10~ s 
0.114149934878319- 10~ 8 



Table 5: The results of the |M| 2 and the sum of all dipole terms are shown. The two implementa- 
tions agree at least to 14 digits for the matrix elements squared and at least to 12 digits for the sum 
of the subtraction terms. The phase space point used in the last entry corresponds to the first one 
in . /Iib/check/check5_uux_w+w-gg/inputm_uux_w+w-gg . h. 
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C-2 


C-l 


CO 


gg -> tig 


AutoDipole 

Ref. oa 


2.49467966948003 • 10~ 4 
2.49467966948004(1)- 10~ 4 


3.68989776683705 • 10~ 4 
3.68989776683706(1) • 10~ 4 


-4.05387364353899 • 10~ 4 
-4.05387364353900(1) • 10~ 4 


uu — > ttg 


AutoDipole 
Ref. EH 


1.38499897972387 • 10" 3 
1.38499897972387(0) -lO" 5 


2.88738914389178 • 10" 5 
2.88738914389179(1) -lO" 5 


-1.56576469322102 • 10" 5 
-1.56576469322102(0) • 10" 5 


ug —* ttu 


AutoDipole 
Ref. m 


3.84580760674706 • 10" 6 
3.84580760674706(0) -10- 6 


7.73777480040817 • 10" 6 
7.73777480040817(1) • 10~ 6 


-5.19929995897616- 10" 6 
-5.19929995897616(1) • 10~ 6 


gu — > ttu 


AutoDipole 
Ref. [[31 


6.22738241305372 -10 _i 
6.22738241305372(0) • 10" 5 


6.81530745255038 ■ 10~ 5 
6.81530745255037(1) • 10" 5 


-1.52377227863896 • 10~ 4 
-1.52377227863896(0) • 10~ 4 



6: The coefficients for the color and spin averaged results for the I-operator. At least 14 digits agreements are obtained. 



